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ABSTRACT 

We study the thermal evolution of primordial star-forming gas clouds using 
three-dimensional cosmological simulations. We critically examine how assump- 
tions and approximations made in calculating radiative cooling rates affect the 
dynamics of the collapsing gas clouds. We consider two important molecular hy- 
drogen cooling processes that operate in a dense primordial gas; H 2 line cooling 
and continuum cooling by H2 collision-induced emission. To calculate the opti- 
cally thick cooling rates, we follow the Sobolev method for the former, whereas 
we perform ray-tracing for the latter. We also run the same set of simulations 
using simplified fitting functions for the net cooling rates. We compare the sim- 
ulation results in detail. We show that the time- and direction-dependence of 
hydrodynamic quantities such as gas temperature and local velocity gradients 
significantly affects the optically thick cooling rates. Gravitational collapse of 
the cloud core is accelerated when the cooling rates are calculated by using the 
fitting functions. The structure and evolution of the central pre-stellar disk are 
also affected. We conclude that physically motivated implementations of radia- 
tive transfer are necessary to follow accurately the thermal and chemical evolution 
of a primordial gas to high densities. 
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Introduction 



The first stars fundamentally transform the early universe by emitting the first light and 
also by synthesizing and dispersing the first heavy elements. They initiat e cosmic reioniza- 
tion, a nd set the scene for the subsequent formation of the first galaxies (see lBromm &: Yoshida 
( 1201 ll ) for a review). Understanding how and when the first stars were formed is one of the 
important goals of modern astronomy. 

There has been a significant progress over the past years in the theoretical studies on 
the first stars. With the currently available computer power, on e can perform an ab initio 



simulatio n of the first star fo rmation in substantial details (see iBromm et al.l ( 120091 ) for a 



review). lYoshida et al] ( 20081) studied t he formation of a primordial protostar in a proper 
cosmological context. iTurk et al.l ( 120091 ) showed that a primordial gas cloud c an fragment 



i nto m ultipl e clumps by the a ction of the rotation and turbulence in the cloud. I Clark et al. 



( 120111 ) and iGreif et al.l ( 1201 ll ) used sink-particle techniques to follow the evolution of the 
accretion disk around a primordial protostar for hundreds years. It was shown that the 
circumstellar disk around a protostar becomes gravitationally unstable to form multiple 
protostars. Unfortunately, these simulations could not be run long enough to obtain the 
solid prediction for multiplicity and for the characteristic mass of the first stars. Including 
the so -called proto- stellar feedback effects is important to determine the mass of the first 
stars ( iMcKee fc Tanll2008l ; IStacy et al.l 120121 ) . Recently, iHosokawa et al.l (120111 ) performed 
radiation-hydrodynamical calculations to show that the self-regulating proto-stellar feedback 
halts the growth of a primordial protostar when its mass is several tens of solar-masses. It 
has become possible to follow the entire evolution of a primordial protostar to the main- 
sequence and thus to discuss rigorously important issues such as the characteristic mass of 
the first stars. 

There still remain a few uncertainties in these theoretical studies on primordial star 
formation. For example, some of t he important chemical reaction rates in a pre-stellar gas 
are not known to good accuracies (IGloverl 120081 ). Most importantly, the three-body hydro- 
gen molecule formation rate is poorly determined, which lea yes substantial u ncertainties in 
the thermal evolution of a pre-stellar gas at high densities ( ITurk et al.l 120 111 ). Calculating 
radiative cooling rates at such high densities, where the gas is optically thick, is essentially 
a radiative transfer problem. Pre yious studies adopt two m ethods to solve this problem. 
Fitting functions are proposed by iRipamonti fc Abel which describes the net cool- 

ing rate as a function of the local density based on the result of a fully one-dimensional 
(ID) radiative transfer calculation. T he other method adopts the large- velocity gradient ap 



proac h, the so-called Sobolev method (lYoshida et al.ll2006l ; IClark et al.ll2011t IHosokawa et al. 



201 ll ). Similar methods have been used to follow the gas evolution at even higher densities 
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( IRipamonti fc Abell 120041 : lYoshida et al.l 120071 120081 ). It is important to examine whether 



or not the different implementations of radiative cooling calculations produce significantly 
different results. 

In this paper, we critically examine whether or not implementations of the optically- 
thick radiative cooling affect the evolution of a primordial gas cloud. To this end, we run a set 
of three-dimensional (3D) cosmological hydrodynamical simulations. We explicitly compare 
the results obtained from the simulations and study the structure of the gas cloud in detail. 

The rest of the paper is organized as follows. In Section 2, we describe the main physical 
processes in a primordial gas cloud evolution. There, we also describe the calculation methods 
of the net cooling rates. The simulation settings and computational methods are given in 
Section 3. Section 4 shows results of the cosmological simulations. We summarize the results 
and give concluding remarks in Section 5. 



2. Thermal evolution of a primordial gas cloud 



A gravitationally contracting primordial gas evolves roughly isothermally. Since the 
onset of run-away collapse, the temperature r ises only a factor of about ten whereas the 
density increases over 16 orders of magnitudes (IPalla et al.lll983l ; lOmukai fc Nishilll998l ). In 
a primordial gas, the main coolant is hydrogen molecules (H2). There are two important 
regimes where radiative transfer effects become important. One is at densities 10 8 cm -3 ^ 
nn ~ 10 14 cm -3 , where H 2 line cooling is dominant and the cloud core cools and condenses 
rapidly. Then the cloud core becomes optically-thick t o lines. It is known that th e 
chemo-thermal instability can be triggered in this phase (ISabano fc Yoshiilll977t ISilk! 119831 ). 



The other is at densities 10 14 cm 3 ^ nn ~ 10 17 cm 3 , where the gas is nearly opaque to H2 
line photons but cooling by the collision-induced emission (CIE) becomes efficient. 

In order to compute the net radiative cooling rate in the two regimes, we need to compute 
the opacity for photons in a broad energy range. In principle, the gas opacity depends on a 
number of physical quantities such as the local gas density, temperature, and velocities. 



2.1. H 2 line cooling 

When the gas density exceeds nn ~ 10 8 cm -3 , rapid three-body reactions of H 2 forma- 
tion convert nearly all the hydrogen atoms into hydrogen molecules. Then H 2 line cooling 
becomes highly efficient. As the gas cloud condenses, however, the gas cloud core becomes 
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opaque to H2 lines. Then the H2 line cooling rate is calculated as 

Ari2, thick ^ hv u [ /3 esca p e it ; -A u [ Tl u , 



(1) 



u.I 



where hv u i is the energy difference between the upper level u and the lower level /, /3 eS cape,uZ 
is the escape probability for a photon without absorption, A u i is the Einstein coefficient for 
the spontaneous transition, and n u is the number density of the hydrogen molecule in the 
upper level u. To calculate the escape probability /3 e sca P e, we evaluate the opacity for each 
H2 line as 

Tlu = OtluL , (2) 

where ati u is the absorption coefficient for the transition from I to u level. Calculating the 
characteristic absorption length scale, L, is the remaining task. To this end, we adopt the 
Sobolev method. The Sobolev length along a line of sight is defined as 

T ^thermal 



\dV r /dr\ 



(3) 



where v thermal = (kT/mn) 1 ^ 2 is the typical thermal velocity of the hydrogen molecules and 
V r is the fluid velocity in the direction. The escape probability for a spherical cloud is given 
by 

1 — exp(— r) 



escape 



T 



(4) 



We compute the Sobolev length and the escape probability in three orthogonal directions 
and use the average value of them to calculate the net cooling rate. 

The above method requires the evaluation of t he optical depths f o r a fe w hundred 
molecular lines, and thus is computationally costly. iRipamonti fc Abel! (120041 ) propose a 
fitting formula for the optically-thick cooling rate based on the result of ID radiative transfer 
simulations. It is given by 



Ah 2 , thick — A H2 ,thin x m i n 



8 x 10 9 cm" 3 



-0.45 



(5) 



Note that the simple function depends only on the local gas density. We use the formula as 
an alternative method to compute the net cooling rate. 



2.2. Collision-induced emission cooling 



At densities greater than ~ 10 14 cm 3 , hydrogen molecules collide so frequently that 
each collision pair temporarily induces an electric dipole and either molecule makes an energy 
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transition by emitting a photon. This process is known as collision- induced emission (CIE), 
which acts as an efficient radiative cooling process at such high densities. The resulting 
emission s pectrum appears essent i ally as continuum rad iatio n. We use the c ollision cross- 
sections of IJorgensen et al.l ( 120001 ) . iBorysow et al.l (j200l[ ). and iBorysowl (120021 ). 



At even higher densities riu > 10 16 cm" 3 , the gas cloud beco mes opaque to the continuum 
emission. We use the Planck opacity table (ILenzuni et al.lll99ll ) to calculate the gas opacity 
and the resulting radiative cooling rate. By noting that the net energy transfer rate should 
scale as A oc 1/(1 + r) for small r, whereas it scales as A oc 1/r 2 for large r, we assume a 
simple form of double power-law 

f = [l + r(n,T)][l + (r(n,r)/10)] (6) 

as the "efficiency" factor of the CIE cooling. Although the particular functional form is 
somewhat ad hoc, it reproduces well the net cooling rate that is obtained fro m more detailed 
radiative transfer calculations (jOmukai Nishil Il998t lYoshida et al.ll2008l ). Note that we 
explicitly write as r = r(n, T) to express that the local opacity depends on the gas density 
and temperature. In practice, we compute the opacity along a line-of-sight by the integral 



k(p, T)p dl 



(7) 



where k is the absorption coefficient and p is the local gas density. In order to take the cloud 
geometry and structure into account, we take the mean of the efficiency factors in three 
orthogonal directions (see equation [6]), 

fx~\~fy~\~ fz /„\ 



/meai 

O 

Then the CIE cooling rate is calculated as 

AciE ,thick AciE,thin x /mean 



(9) 



A simple continuum opacity model is also proposed by lRipamonti fc Abel! (120041 ) . which 
is given by a function of the local gas density nn 2 as 



AciE,thick — AciE.thin x mm 



-TCIE 



where 



TCIE 



n U2 



7 x 10 15 cm- 3 



TCIE 



2.8 



(10) 



We compare the runs with the above two methods of computing the CIE cooling rate, 
(as equation [9] and fTU]). 
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Numerical Simulations 



We use the p arallel iV-body/Smoothed Particle Hydrodynamics (SPH) solver GADGET- 
2 ( Springelll2005l ) in its version suitably adopted for the primordial star formation. We solve 
chemical rate equations for fourteen species of primordial species (e~, H, H + , He, He + , 
He ++ , H 2 , H j~, H~, D, D + , HD, HD + , HD~). The reactions and rates are summarized in 



Yoshida et all (120061 . 120071 ). We adopt the A-Cold Dark Matter (A-C DM) cosmology. Th e 
cosmological parameters are based on the seven-year WMAP results (ILarson et al.ll201ll ); 
Vl A = 0.734, tt m = 0.236, £l b = 0.0449, and H = 71.0 km s" 1 Mpc" 1 . The normalization 
of the power spectrum is set to be erg = 2.0 such that structure forms early in the small 
simulation volume. All simulations are initialized at z ini = 99. 

We use the zoom-in re-simulation techniques to achieve a large dynamic range. First, 
we run a cosmological simulation to locate a star-forming halo which is later re-simulated 
with a higher resolution. The comoving box size of this parent simulation is 50 kpc h _1 on 



a side. In the zoomed region, the initial particle masses are m gas 
mdm = 1.01 x 10~ 2 Mq, respectively. 



2.04 x 10" 3 M and 



We run each of the zoomed simulations first to the point where the central density 
reaches n cen = 10 8 cm" 3 . Then, to cal c ulate the gas collapse further, we use the particle 
split method of iKitsionas fc Whitworthl (120021 ) to achieve a higher mass resolution. We do 
this refinement progressively such that a local Jeans length is always resolved by 10 times 
the local SPH smoothing length. By using this technique, we achieve a mass resolution of 
4 x 10~ 6 M Q at the last output time. We stop the simulations when the central 



m 



Wis 



density reaches n c 



10 17 cm" 



4. Results 

To examine the properties of star-forming halos, we choose two realizations as character- 
istic cases, Run A and B. Run A forms a disk-like structure inside the collapsing gas cloud, 
whereas Run B shows elongated structure and eventually develops S-shape arms (see the top 
panels in Figures [Hand [2]). The bottom panels show, in both cases, the central sub-parsec re- 
gion is significantly flattened. The non-spherical structures likely yield direction-dependence 
of both line and continuum opacities. 

In the following subsections, we first compare the opacity models for H 2 line cooling. 
Then we examine differences in the CIE cooling phase. We have already mentione d the 



importance of the three-body molecular hydrogen formation rate (ITurk et al.ll201ll ). We 



examine the effect of varying the three-body reaction rate at the end of this section. 
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4.1. H 2 line cooling 

Figure |3] shows the radial profiles of five physical quantities for Run A when the central 
density is ~ 10 14 cm -3 . At this time, the central part is nearly completely opaque to H2 
line emission. The solid line shows the result from the run with the 3D opacity calculation 
(Sobolev method, equation pQ) and the dotted line is for the run with the fitting opacity 
formula (equation [5]). We see significant differences in the temperature, radial velocity, and 
accretion rate profiles. 

Because the photon escape probability depends on the local velocity gradient (see equa- 
tion [3]), the line cooling rate can be large if the velocity gradient is large in the dense cloud 
core. The radial infall velocity is critically affected by the degree of rotation of the collapsing 
cloud. Star-forming clouds in a cosmological simulation generically have finite initial angular 
momenta and thus they spin up gradually as they collapse gravitationally. The radial infall 
velocity, V r , and the velocity gradient, dV r /dr, in such rotating gas clouds are smaller than 
realized in spherically symmetric collapse. Thus the escape probability in the cosmological 
simulation is smaller than the fitting formula predicts. 

In order to investigate further the difference in the escape probability, we perform a 
spherical collapse simulation with 3D set-up. We follow gravitational collapse of a super- 
critical Bonnor-Ebert sphere having a mass of ~ 1000 M . For this run, we calculate the 
H2 line opacity using the Sobolev method. The results are shown in Figure [3] (dashed line 
in each panel). Clearly, the radial velocity of the run is larger than that of our cosmological 
simulation Run A, which has a substantial degree of rotation. This is the major source of 
the differences in the line escape probability and in the cooling rate, as shown in Figure HI 

The left panel of Figure H] shows the escape probability of H2 line photons as a function of 
the gas density. We show three snapshots for the mean profiles (solid lines) when the central 
density is n cen = 10 10 , 10 12 , and 10 14 cm~ 3 . The escape probability calculated by our 3D 
treatment is smaller than the fitting formula (dotted line) most of the time. The difference 
is as large as a factor of ten at the densest part. The fitting formula over-estimates the net 
cooling rate. This is easily understood by the collapse speed of the spherically symmetric 
calculation. 

Next, let us consider the direction-dependence of the escape probability. In the left 
panel of Figure IU we plot the escape probability in the direction along x— , y—, and z-axes 
(long-dashed, short-dashed, and dot-dashed lines). We configure the coordinate such that 
the ^-direction is aligned to the angular momentum vector of the central cloud core @. The 



Figures Q] and [2] also use the same coordinate. 
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escape probability is large along the z-direction. Thus, physically, line photons preferentially 
escape in perpendicular directions to the flattened cloud core. 

It is important to note that the difference (or mis-estimate) in the cooling rate critically 
affects the collapse dynamics. The right panel of Figure H] shows the time evolution of the 
central density since the central density reaches n cen = 10 8 cm -3 . In the run with the fitting 
formula, the cloud core collapses earlier by ~ 20, 000 years, i.e., the gravitational collapse is 
accelerated because of the "efficient" cooling. 

Figures [5] and [6] show the same results for Run B, which has a spiral structure. The 
overall evolutionary trend is quite similar to Run A. Clearly, the importance of radiative 
transfer effects is not particular to the configuration of Run A. The differences in the radial 
profiles of Run B are understood similarly to Run A as explained in the present section. 
We conclude that the multi-dimensional treatment for the radiative cooling is important to 
follow the thermal evolution and the gravitational collapse accurately. 

4.2. H 2 collision-induced emission cooling 

Next, we discuss the thermal evolution through the phase where CIE cooling is impor- 
tant. In this subsection, we primarily discuss the results of Run B. The radial profiles when 
the central density is n ccn ~ 10 17 cm -3 are shown in Figure [7J We also plot the normalized 
cooling rate, the efficiency factor / = AciE,thick/AciE,thin and the time evolution of the central 
density in Figure [HJ 

The left panel of Figure El shows similar features to the case of H 2 line opacity discussed 
in the previous section. We a significant direction-dependence of the escape probability. 
The fitting formula over-estimates the net cooling rate compared with the run with our 3D 
opacity treatment. The difference in the net cooling rate becomes as large as a factor of 5 
at ^ccn ~ 5 x 10 15 cm" 3 . The gas collapses faster with the fitting formula for the opacity 
calculation. However, the difference in the resulting collapse time is not very large (see 
the right panel of Figure [8]). The time difference is only about 1 year at the end of the 
calculations. 

It is worth pursuing the reason why the collapse proceeds similarly in the two cases 
despite the large difference in the net cooling rate in the relevant regime nn ~ 10 15 — 
10 16 cm -3 . To this end, we perform two additional calculations for Run B. One is run with 
an artificially increased CIE cooling rate as Acie ->5x Acie whereas the other is run with 
the reduced efficiency factor / mcan -4 0.2 x /mean- The other configurations are identical 
to Run B, and the opacity calculations are also done using 3D ray-tracing. Therefore, we 
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expect that the two runs clarify the impact of enhanced/reduced CIE cooling. The results are 
shown in Figure |9j In the case with the increased cooling rate, the collapsing gas evolves on 
a lower temperature track (the long-dashed line in Figure ED, as is naively expected from the 
enhanced cooling rate. Clearly, the cooling rate itself directly affects the gas temperature in 
this regime. On the other hand, the effect of decreasing the efficiency factor appears relatively 
small (the short-dashed line). The cloud evolves on a slightly higher temperature track, but 
the difference is significant only in a narrow range of density, 10 15 < nn < 10 16 cm" 3 , where 
the gas is collapsing rapidly. The free-fall time there is estimated to be ~ 1 year which 
is comparative to the cooling time. The core condenses quickly and becomes optically thick 
to the continuum photons after the central density reaches ~ 10 16 cm -3 . In the left panel 
of Figure we also plot the time evolution of the cooling efficiency at the central part 
(double-dotted line). The evolution looks similar to the fitting function at the high density. 

In summary, the accuracy of the continuum opacity calculation causes only minor effect 
on the thermal evolution. The difference in / = AciE,thick/AciE,thin is large between the 
methods, but the resulting evolution is not sensitive to the details of the methods. It should 
be noted, however, that we have considered only a particular case in which the central core 
undergoes rapid run-away collapse. In other circumstances, for example in an accretion disk 
around a protostar where the density evolution is much slower than in the collapsing gas 
core that we have studied, accurate calculations of the optically-thick cooling may be more 
important. 



4.3. Three-body H 2 formation 

It is well known that there is a large uncertainty in the reaction rate of the three-body 
H2 formation 

H + H + H^H 2 + H. (12) 

Because this is the dominant reaction to form hydrogen molecules at high densities (tih > 
10 8 cm -3 ), an accurate reaction r ate is needed to determine the chemical and thermal evo- 



lution of a primordial gas cloud. iGloverl ( 120081 ) summarize various rate coefficie nts used in 



the lit erature, which differ by a factor of 30 at the relevant temperature range. iTurk et al. 



( 1201 ll ) perform a set of hydrodynamical simulations to directly study the overall effect caused 
by the uncertainties in the reaction rate. They conclude that, while the difference between 
different realizations (gas cloud samples) is larger than that caused by the uncertainty of 
the three-body rate coefficient, the morphology and the collapse time of a gas cloud depend 
strongly on the reaction rate. 



In this section, we revisit the issue because the density range where the uncertainties 
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are relevant is coincident with the range where the radiative transfer treatment is important, 
as studied in the previous sections. We are able to compare the overall differences caused 
by the uncertainties of the reaction rates with the differences caused by radiat ive transfer 



treatm ents. So far in the present paper, we have adopted the reaction rate from iPalla et al. 



( 1l983l ) (Medium case in Figure [TO]) . We run the same sim ulation as R un B but w i th the 
different three-body reaction rates of iFlower fc Harris! (120071 ) (High) and lAbel et al . (2002) 
(Low) , which are the largest and smallest reaction rates among those compiled by I Glover 
(120081 ). The respective rates are summarized in Figure [TU1 and Table [TJ 

Figure [11] shows the simulation results with the three reaction rates. The left panel 
shows the H2 fraction as a function of gas density. The density at which the cloud becomes 
fully- molec u lar d iffers more than a factor of 10. This is consistent with the conclusion of 



Turk et al.l (1201 ll ) . Because the molecular fraction largely determines the H 2 line cooling 
rate, the resulting collapse time of the cloud differs by At ~ 30, 000 years, depending on the 
choice of the reaction rate. It is interesting that the time difference of the cloud collapse, At, 
is comparable to the difference caused by the choice of the opacity calculation method. The 
difference in the molecular fraction becomes large at > 10 9 cm -3 , where the calculation of 
the H 2 line opacity is important (see Figure E]). Therefore, we argue that using an accurate 
radiative transfer method is as important as using an accurate reaction rate. 



5. Discussion and Conclusion 

Radiative cooling by hydrogen molecules governs the thermal evolution of a primordial 
star-forming gas cloud. At high densities, the cloud core becomes optically thick to both 
the rot-vibrational lines and collision-induced continuum. It is necessary to estimate the gas 
opacities in multiple directions in order to calculate the radiative cooling rates accurately. We 
have explicitly compared the results from several sets of simulations with different manners 
to calculate the optically-thick cooling rates. 

When non-spherical structures develop in the central region of the cloud, photons do not 
escape isotropically from the dense part. Our simulations show that the gas cloud spins up as 
it contracts, and forms a flattened disk at the center. Then the photon escape probability not 
only varies with time but is also direction-dependent. Utilizing an "isotropic" fitting formula 
that is derived from spherically symmetric calculations over-estimates the net cooling rate 
and causes the cloud core to collapse fast (see Figures H] [6] and [8]). With our 3D opacity 
calculation, the photon escape fraction from the cloud core is always smaller than given by 
the fitting formula. The resulting cooling rate differs by a factor of a few to 10, depending on 
the exact density, temperature, and velocity structure. There is also a directional effect of 
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the radiative transfer. In perpendicular directions to the faces of the disk-like cloud, photons 
can easily escape. 

Details of the implementation of the optically-thick radiative cooling affects the thermal 
and dynamical evolution of the cloud. Figure [12] shows the density distribution of two 
characteristic runs. We use the snapshots at 10 years after the central density reaches 
n cen = 10 18 cm -3 . The left panel shows the case with our 3D radiative transfer, whereas 
the right panel is for the case with the fitting formulae q The latter case appears rounder 
and more concentrated than former, which has a spiral structure. Clearly, further detailed 
studies on radiative transfer effects are needed, particularly on the long-term evolution of 
the proto-stellar disk. The exact structure of the proto-stella r disk likely affec ts the disk 
evolution, mass accretion rate and fragmentation of the cloud (IGreif et al.ll2012[ ). 



A similar comparat ive study of radiative cooling implementations has been done by 
Wilkins fc Clarkd (120121 ) in the context of the present-day star formation . Interestingly, an 
oppos ite trend is found in the case of "polytropic" cooling examined by IWilkins fc Clarke 
( 120121 ) . which is based on locally estimated opacities. They show that the polytropic cooling 
performs well only in spherically symmetric cases. The polytropic method over-estimates the 
column density, and hence underestimates the radiative cooling rate in non-spherical cases. 
In order to calculate radiative cooling rates in 3D simulations, it is important to take the 
direction-dependence of the photon diffusion into account, similarly to what we conclude in 
the present paper. 

It is clearly advantageous to use a computational method that is fast and robust. Our 
method of the H2 line transfer utilizes local velocity gradients that come with essentially 
no additional cost, because the velocity gradients are already computed and used in the 
other parts in our smoothed-particle hydrodynamics code. For continuum photons, we need 



to compute the co lumn density a 



ong six (or more) directions using a costly projection 



method devised by lYoshida et al.l (120071 ). In fact, the continuum opacity calculation is one 



of the most time consuming part in our 3D simulations. Nevertheless, we argue that it is 
necessary to properly take the direction-dependence into account in order to calculate the 
optically-thick radiative cooling rate accurately. We also note that our method is based on 
the so-called escape probability method, which itself is an approximation. Essentially, we 
assume only the densest part emits continuum photons. It is desirable to implement fully 
three-dimensional radiative transfer, by employing advanced methods such as flux-limited 



2 These calculations are performed with a low resolution such that the central part does not collapse to 
much greater than 7jh = 10 18 cm~ 3 . We keep the low mass resolution deliberately in order to follow the disk 
evolution over 10 years. 



diffusion or Ml-closure (e.g., IWhitehouse fc Batd (12004 ); 
long-term evolution of a primordial proto-stellar system. 



Levemore 



([19841)) to follow the 



Hydrodynamical simulations with radiative cooling are commonly used for the study of 
the primordial star formation. In such simulations, it is important to use accurate methods to 
calculate radiative cooling rates. For example, whether or not a proto-stellar disk fragments 
is determined by the thermal and gravitational instability of the circumstellar gas. The 
disk fragmentation is an important issue which is thought to determine the multiplicity and 
possibly the characteristic mass of primordial stars. Our study clarifies the importance of 
mult i- dimensional radiative processes in a primordial star-forming cloud. 

We thank Hideyuki Umeda for stimulating discussions and Ikko Shimizu for the tech- 
nical support. The work is supported in part by the Grants- in- Aid for Young Scientists 
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YITP in Kyoto University and T2K-Tsukuba System in Center for Computational Sciences, 
University of Tsukuba. 
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Fig. 1. — Projected density distribution for Run A, which shows a flattened disk-like struc- 
ture. Face-on views (top panels) and edge-on views (bottom panels) in a volume of, from 
left to right, 1, 10 _1 , 10~ 2 , and 1CT 3 pc ~ 200 AU on a side, respectively. 
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Fig. 2. — As for Figure [U but for Run B. 
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Fig. 3. — Radial profiles of various quantities for Run A. Five panels show the enclosed 
mass, temperature, radial velocity, gas number density, and gas mass accretion rate, in the 
clockwise order starting from the top-left. The solid line shows the result with the 3D H2 line 
opacity calculation (Sobolev method), the dotted line is for the run with the fitting function 
and the dashed line is for a spherical collapse calculation with 3D set-up. 
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Fig. 4. — Left panel: The effective escape probability, A t hi C k/Athin, for H 2 line cooling for 
Run A. The long-dashed, short-dashed, and dot-dashed lines show the three orthogonal 
components (x, y, and z) when the central density is n cen = 10 14 cm -3 . The solid lines 
show the evolution of the direction-averaged mean value at n cen = 10 10 , 10 12 , and 10 14 cm -3 . 
The dotted line shows the fitting function given by equation flH]). The double-dotted line 
indicates the result of a spherical collapse calculation with 3D treatments. Right panel: The 
time evolution of the central density. The horizontal axis is the elapsed time after the central 
density reaches ~ 10 8 cm -3 . 
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Fig. 5. — As for Figure [3] but for Run B. 
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Fig. 6. — As for Figure H] but for Run B. 
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Fig. 7. — Radial profiles of various quantities for Run B in the high density regime where 
CIE cooling is important. Five panels show the enclosed mass, temperature, radial velocity, 
gas number density, and mass accretion rate, in the clockwise order starting from the top- 
left. The solid line is for the result with the 3D H 2 CIE opacity calculation, the dotted line is 
for the run with the fitting function and the dashed line is for a spherical collapse calculation 
with 3D treatments. 
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Fig. 8. — Left panel: The net cooling rate by H2 CIE emission for Run B. The long-dashed, 
short-dashed, and dot-dashed lines show the three orthogonal components (x, y, and z) when 
the central density is n cen = 10 17 cm -3 . The solid lines show the mean values of them at 
n ccn = 10 15 , 10 16 , and 10 17 cm -3 . The double-dotted line is the time-evolution of the central 
part. The dotted line shows the fitting function given by equation f fTOj) . Right panel: The 
time evolution of the central density. The horizontal axis is the elapsed time after the central 
density reaches ~ 10 14 cm -3 . 




Fig. 9. — Thermal evolution in the temperature - density plane for Run B. We plot the 
simulation results with our 3D radiative transfer treatment (solid), with the fitting opacity 
function (dotted), with a increased CIE cooling rate (long-dashed), and with a reduced CIE 
escape probability (short-dashed). The thick lines show the final output time and the the 
thin lines are for the results at earlier phases. 
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Fig. 10. — Three-bo dy formation reactio n rate used in the literature. Th e solid line 



(High) is the rate from lFlower fc Harris! (12007 ), the d ashed line (Medium) is from lPalla et al. 



( 119831 ) , and the dotted line is from lAbel et al. 

m 



( 120021 ). The rates are explicitly given in Table 
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Fig. 11. — Left panel: Molecular fraction as a function of the gas density. Three lines 
show the runs with three reaction rates (see Figure [10] and Tabled]). Right panel: The time 
evolution of the central density. The horizontal axis is the elapsed time after the central 
density reaches ~ 10 7 cm" 3 . 
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Fig. 12. — Projected density distributions for the same realization but with different opacity 
calculations at 10 years after the central density reaches 10 18 cm -3 . The plotted region is 
200 AU on a side. The left panel is for the result with our 3D opacity calculation, whereas 
the right panel is for the run with the fitting opacity for line emission. 
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Table 1. Three-body H2 formation rate 



Model 


Rate Coefficient (cm 6 s 1 ) 




Reference 


High 


1.44 x 1CT 26 T -1 - 54 




1 


Medium 


5.5 x 10- 29 T- L0 




2 


Low 


1.14 x IO- 31 T-°- 38 


(T < 300 K) 


3 




3.9 x 1(T 30 T- 10 


(T > 300 K) 


3 


Note. - 


- T is the gas temperature in 


Kelvin. 
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